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Abstract. We present time-domain Monte Carlo simulations of radio emission from cosmic ray air showers 
in the scheme of coherent geosynchrotron radiation. Our model takes into account the important air shower 
characteristics such as the lateral and longitudinal particle distributions, the particle track length and energy 
distributions, a realistic magnetic field geometry and the shower evolution as a whole. The Monte Carlo approach 
allows us to retain the full polarisation information and to carry out the calculations without the need for any 
far-field approximations. We demonstrate the strategies developed to tackle the computational effort associated 
with the simulation of a huge number of particles for a great number of observer bins and illustrate the robustness 
and accuracy of these techniques. We predict the emission pattern, the radial and the spectral dependence of the 
radiation from a proto typical 10 17 eV vertical air shower and find good agreement with our analytical results 
jHueee fc Falckel 120031) and the available historical data. Track-length effects in combination with magnetic field 
effects surprisingly wash out any significant asymmetry in the total field strength emission pattern in spite of the 
magnetic field geometry. While statistics of total field strengths alone can therefore not prove the geomagnetic 
origin, the predicted high degree of polarisation in the direction perpendicular to the shower and magnetic field 
axes allows a direct test of the geomagnetic emission mechanism with polarisation-sensitive experiments such as 
LOPES. Our code provides a robust, yet flexible basis for detailed studies of the dependence of the radio emission 
on specific shower parameters and for the inclusion of additional radiation mechanism in the future. 

Key words, acceleration of particles - elementary particles - polarization - radiation mechanisms: non-thermal - 
methods: numerical 

classical particle detecto r techniques at moderate cost 
l|Falcke fc Gorhamll2003|) . Out of this interest, the LOPES 
project was initiated as a joint effort to firmly establish the 
feasibility of radio measurements o f extensive air showers 
with a LOFAR prototype station ijHorneffer et all l2003|) 
operating in conjunc tion with the KASCADE experiment 
(Ant oni et al.U2003(l as well as to gain a solid theoretical 
understanding of the emission mechanisms. 

In iHuege fc Falcke] l)2003j) we presented calculations 
of radio emission from extensive air showers within the 
scheme of "coherent geosync h rotron radiation" first pro- 
posed by iFalcke fc Gorhaml l|2003|h These calculations 
were based on an analytic approach and were specifically 
aimed at gaining a solid understanding of the coherence 
effects that shape the radiation emitted by an air shower. 
Building on this foundation, we now continue to develop 
and enhance our model with elaborate Monte Carlo (MC) 
simulations. The MC technique allows us to infer the emis- 
sion characteristics with much higher precision by taking 



1. Introduction 

Extensive air showers initiated by high-energy cosmic rays 
have been known for almost 40 years to produce strongly 
pulsed radio e mission at frequenc ies from a few to a few 
hu ndred MHz (Jellev et aT1ll965t for an excellent review 
sec lAllanl Il97ll) ~ The quality of the experimental data 
as well as the theoretical understanding of the emission 
mechanism, however, stagnated on a rather basic level 
when the research in this field almost ceased completely 
in the late 1970ies. 

Lately, the advent of new fully digital radio inter- 
ferometers such as LOFAR 1 sparked renewed interest 
in the topic. The radio technique potentially offers a 
cost-effective additional method for the study of exten- 
sive air showers which would very much complement the 

Send offprint requests to: T. Huege, 
email: thuege@mpifr-bonn.mpg.de 
1 http://www.lofar.org 
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into account more realistic and complex shower properties 
and applying fewer approximations than in the analytic 
calculations. At the same time, it provides an independent 
means to verify our previous calculations due to the to- 
tally different technique employed. The app roach we take 
is sim ilar to the MC simulations done by ISuprun et alJ 
l)2003(l . yet our simulation is developed to a much higher 
level of complexity. 

The layout of this paper is as follows: In sections 2 
and 3 we motivate and explain the application of the MC 
technique and provide details about its implementation. 
In section 4, we explain the "intelligent" concepts that we 
have explored to make simulations with a high number of 
particles and observer bins feasible on standard personal 
computer hardware. After a short description of the raw 
output of our program and the associated data reduction 
in section 5, we demonstrate the robustness and consis- 
tency of the implemented algorithms in detail in section 
6. Similar to our earlier theoretical calculations, we then 
first concentrate on the emission from a single "slice" of 
particles in the air shower and compare the results with 
the analytical results in section 7. Next, we perform the in- 
tegration over the shower evolution as a whole and present 
the results in comparison with our analytical work and his- 
toric data in section 8 before concluding with a discussion 
of the results and our conclusions in sections 9 and 10, 
respectively. 



2. The Monte Carlo approach 

After having established the general dependencies of the 
radio emission on a n umber of air shower properties in 
iHuege fc Falcke] (j2003), a continuation of the calculations 
with MC techniques offers a number of advantages. 

2.1. Motivation and objectives 

First, MC techniques allow an independent verification 
of the analytic calculations by adopting the exact same 
shower properties (e.g., distributions of particles in space 
and energy, choice of magnetic field and shower geome- 
try), but doing the calculations in a completely different 
way. In particular, the MC calculations are carried out by 
summing up the individual particles' pulses in the time- 
domain, whereas the analytic calculations were done in 
the frequency domain. 

Second, it is relatively easy to include highly complex 
(and thus realistic) shower characteristics in the MC sim- 
ulations, thereby increasing the model precision signifi- 
cantly over the previous results. Additionally, it is not 
necessary to make approximations such as adopting a far- 
field limit, which reduces the accuracy of the analytic cal- 
culations. 

Overall, MC simulations therefore constitute the logi- 
cal next step in the development of our model. 



2.2. General approach 

The general idea of a MC simulation of radio emission 
from cosmic ray air showers is simple: 

1. model the radiation emitted by an individual particle 
as precisely as possible, 

2. distribute particles randomly in a simulated shower 
according to the desired shower characteristics (e.g., 
spatial and energy distributions), 

3. superpose the radiation received from the shower par- 
ticles at the given observing positions, taking into ac- 
count retardation effects. 

In fact, the "microphysics" of an individual particle's 
geosynchrotron emission (step 1) can be described analyt- 
ically without the need of any approximations as long as 
one does not simulate particle interactions explicitly but 
only takes them into account via statistically distributed 
track lengths. The strategy for the implementation of steps 
2 and 3 then is as follows: 

— generate shower particles according to the desired dis- 
tributions, 

— for each particle: 

— for each observing position (ground-bin): 

• establish an adequate sampling of the particle 
trajectory, 

• calculate and retard the emission contributions 
emanating from the sampled points on the tra- 
jectory, building up the electric field time-series 
that the observer sees, 

• incorporate the contributions into the ground- 
bin's pre-existing time-series data. 

This strategy is fairly simple. What makes the problem 
difficult is the huge computational effort of a simulation 
with a large number of particles and g round-bins, as was 
already discussed by iDova et all l|l999j) who considered a 
similar approach. A number of intelligent concepts has to 
be applied to actually develop a working simulation out of 
this simple "recipe". We discuss these concepts in depth 
in section 21 



3. Implementation details 

In this section we describe a number of relevant imple- 
mentation details of our MC code. 

3.1. Technical information 

The program has been developed under Linux using gcc 
2.9.5, gcc 3.3.1, and the Intel C++ compiler version 8.0.5 
which produces much faster code. No third-party libraries 
were used in order to maximise portability and minimise 
dependency on external factors. The program source-code 
will be made available at a later time. 



T. Huege and H. Falcke: Radio emission from cosmic ray air showers 



3 



3.2. Particle creation and propagation 

The particles in the shower are created with random prop- 
erties distributed according to analytic parametrisations. 
If not explicitly stated otherwise, the param etrisations 
are chosen exactly as in lrlnesre fc Falcke] (|2003fl . to which 
we refer the reader for detailed information. While these 
parametrisations are admittedly crude and do not take 
into account some air shower properties such as a realis- 
tic particle pitch-angle distribution, this approach allows 
a direct comparison of the MC results with our analytic 
calculations. A more realistic modelling of the air shower 
will be achieved once our code is interfaced with the air 
shower simulation code CORSIKA. The particles are cre- 
ated with the following properties chosen randomly: 

— the shower age at which a par ticle is created (longitudi- 
nal development according to lOreisenl l|l960l) function), 
which directly yields 

— the position along the shower axis 

— the creation time 

— the lateral shift from the shower axis (NKG- 
distrib utio ri dating back to iKamata fc Nishimural 
<|l958l) and lOreisenNl960l^ 

— the longitudinal shift along the shower axis as function 
of the lateral shift (asymmetrical T-PDF) 

— the azimuth angle for the lateral shift (isotropic) 

— the particle gamma factor (broken power-law distribu- 
tion or fixed 7 = 60) 

— the particle track length (exponential probability dis- 
tribution or fixed A = 40 g cm~ 2 ) 

To take into account the pair- wise creation of particles, 
one electron and one positron are always generated with 
the same properties. At the moment, no random spread in 
the particle pitch-angle is introduced, i.e., the initial par- 
ticle momenta radially point away from a spherical surface 
with 2,300 m radius, exactly as in the analytical c alcula - 
tions, motivated by the data from lAgnetta et all l|l997j) . 
The fact that the initial velocity direction is shared by 
both the electron and positron introduces only minor error 
as the transverse momentum arising from the pair produc- 
tion is minimal — a fact that is intuitively illustrated by 
the still very dense core of the lateral distribution function 
even after many generations of particle creation. 

As each of the particles is created at a specific position 
at a given time with a given initial velocity, the trajectory 
r(t) resulting from the deflection in the geomagnetic field 
is a well-defined helix which can easily be described ana- 
lytically by equations , (|2fJ|> and as derived in the 
appendix. 



3.3. Calculating and collecting contributions 

Once the trajectory (and its time-derivatives) are known 
analytically, the radiation an observ er at position x re- 
ceives at time t can be calculated (cf. IJacksonlll975l equa- 



tion 14.14) immediately by 



E{x,t) 



n-/3 



7 2 (1 -/3-n) 3 i? 2 

n x {(n - (3) x /3} 
(1 -(3-nfR 



(1) 



where e denotes the particle charge, /3(t) = v(t)/c is di- 
rectly given by the particle velocity, R(t) refers to the vec- 
tor between particle and observer position, R(t) = \R{t)\ 
and n{t) = R{t)/R{t) denotes the line-of-sight direction 
between particle and observer. 

The index "ret" points out that the quantities in the 
brackets have to be evaluated at the retarded time 



trot = t - i?(i rct )/c 



(2) 



rather than at the time t in order to accommodate the fi- 
nite light-travel time. This "recursive" retardation relation 
imposes significant problems for an analytical calculation 
in the time-domain. In case of a MC simulation, on the 
other hand, it is absolutely straight forward to take the re- 
tardation into account by simply delaying the emitted sig- 
nal appropriately before collecting it in the ground-bins. 

The first term in equation constitutes the "static" 
term that falls off with R~ 2 in field strength. It is 
usually neglected, an d was not taken into account in 
iHuege fc Falcke) (2003) either. While its contributions are 
indeed negligible, we still include it in our MC simulation 
as to not make any unnecessary approximations. 

The second term is the usual "radiation" term which 
drops as R^ 1 in field strength and therefore dominates 
very quickly over the static term as one goes to higher 
distances. In analytic approaches, it is usually neces- 
sary to apply approximations such as the Fraunhofer- 
approximation for the far-field limit to this term, which 
naturally limits the precision of the results. Again, for a 
MC simulation, it is not necessary to make any approxi- 
mation for the radiation formula. 

Beaming effects are naturally taken into account in 
this formula through the (1 — (3 ■ n) 3 terms in the denom- 
inator that lead to very high field strengths for particle 
velocities close to c and small angles to the line-of-sight. 
As soon as one takes into account the refractive index of 
a medium rather than vacuum, the denominator actually 
becomes zero at the Cerenkov angle. The arising singular- 
ity, in combination with the modified retardation relation, 
then leads to Cerenkov radiation. The analysis of these ef- 
fects, however, is beyond the scope of this work and will 
be carried out in a later paper. 

For a given particle, the trajectory is then sampled 
in a sufficiently high number of points, and the retarded 
contribution emanating from each of these points is calcu- 
lated for each of the observing ground-bins with the full 
precision of equation (Q. Thus, the retarded electric field 
time-series produced by the particle is inferred for each of 
the ground-bins and can then be incorporated into their 
pre-existing time-series data. The details and subtleties of 
this procedure are explained in sections 14.21 and 14.31 
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Layer 


Height [km] 


at [g cm~ 2 ' 


bi [g cm 2 ] 


d [cm] 


1 


0-4 


-186.56 


1222.66 


994186.38 


2 


4-10 


-94.92 


1144.91 


878153.55 


3 


10 - 40 


0.61 


1305.59 


636143.04 


4 


40 - 100 


0.00 


540.18 


772170.16 



Table 1. Parameters for the parametrisation of the at- 
mospheric layers. 



3.4. Atmosphere model 

We use the US stand ard atmosphe re of 1977 as imple- 
mented in CORSIKA l|UlrichL HooH) . It splits the atmo- 
sphere in a number of layers, in each of which the path 
depth in g cm~ 2 (as counted vertically from the outer edge 
of the atmosphere) is parametrised by an exponential de- 
pendence on height 



X(h) = ai + bi exp 



(3) 



The layer boundaries and corresponding parameters a^, 
bi and Cj are giv en in table [T] Th e Moliere radius tm is 
parametrised as l|Dova et all l2003|) 



tm{X) 



9.6 



(X - oi) 



(4) 



3.5. Random number generation 

An important centre-piece of any MC simulation is the 
random number generator at its heart. We revert to the 
well- known "Mersenne- Twister" rand om-number genera- 
tor of iMatsumoto fc Nishimural l|l998|) in the C++ imple- 
mentation by If possible, the generators for 
non-uniform probability distributions were implemented 
using analytic inversions. If not, we reverted to rejection 
methods. 



4. Intelligent concepts 

As mentioned earlier, a "brute-force" approach as 
sketched in the "recipe" given in section lT^l involves such a 
high computational effort (both regarding CPU time and 
memory) that it simply is not feasible on standard PCs. 
One therefore has to employ a number of intelligent con- 
cepts to minimise the computing effort. The main ideas 
that we have explored to reach this goal are described in 
the following subsections. 



4.1. Cutting off 7 1 -cones 

The radiation pattern emitted by an individual highly rel- 
ativistic parti cle is heavily beamed in the forward direc- 
tion (see, e.g.. I Jacksorll975j) . Most of the emission is radi- 
ated into a cone with opening angle of order ~ 7 _1 - This 
directly leads to a very simple, but effective idea of how 
to cut down on computing time: 




Fig. 1. Cutting off 7 1 -cones: only bins in the "ground- 
trace" of a particle's trajectory are selected for evaluation. 

For an individual particle flying on its given trajectory, 
the 7 _1 emission cone sweeps over a relatively small re- 
gion on the ground. It is only in this ground-region of a 
few times 7 -1 "width", which we call the "ground-trace" 
of the particle trajectory, that receives considerable contri- 
butions of radiation from that particle. Thus, we can select 
only the ground-bins inside the ground-trace of this spe- 
cific particle for evaluation as sketched in figure^ hugely 
cutting down on computing time. 

It turns out, however, that the discrete cutting off after 
a certain angular distance introduces errors in the calcu- 
lation, which are — although only at a percent-level — 
significant when one is interested in the emission strength 
at distances a few hundred metres from the shower core. 
We will discuss the details in section IB~51 

4.2. Smart trajectory-sampling 

In general, there are two approaches of how to calcu- 
late the time-dependence of the emission that a specific 
ground-bin receives from a given particle: 

The first approach is to take the trajectory of a given 
particle, sample it in a sufficient number of points and 
calculate the corresponding contributions for the given 
ground-bin using (Q. If the particle trajectory is sam- 
pled in equidistant time-steps, the retardation effects in- 
volved lead to a non-equidistantly sampled time-series in 
the ground-bin. The hcavily-pcakcd pulse shape is, how- 
ever, automatically sampled with high precision in this 
approach. 

As the time-series data collected by the ground-bin 
have to be gridded eventually, it would be easier if one 
could take another approach: taking the ground-bin's pre- 
defined equidistant time-grid and sampling the particle 
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Fig. 2. Smart trajectory sampling: the particle trajectory 
is sampled densely in the peak-region and only sparsely 
further outside. 
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trajectory in the points corresponding to this grid. This, 
however, is not easily possible because of the "recursive" 
retardation relation (01. One would have to do an iter- 
ative search for the corresponding points on the particle 
trajectory, which would be at least as time-consuming as 
interpolating or binning a non-gridded time-series derived 
with the first approach to the ground-bin's time-grid. At 
the same time, this approach bears the risk of missing 
highly peaked contributions in case of a too widely spaced 
time- grid. 

The first approach, therefore, is the better choice for 
our calculations. It is, however, possible to improve on 
the case of equidistant sampling of the particle trajectory 
by taking advantage of the strong beaming of the emis- 
sion once again. The peaks in the time-series result when 
the denominator (1 — (3 ■ n) 3 in Q gets small as the an- 
gle between (3 and n gets small. It therefore makes sense 
to densely sample the region of small angle to the line- 
of-sight, and increase the distance between the sampled 
points (up to a maximum value) as the angle gets larger. 
This minimises the number of points used while it guar- 
antees high-precision sampling of the pulse shape as can 
be seen in figure El 

4.3. Intelligent gridding strategy 

The individual particles' pulses are very short, of order 
10~ n s for 7 = 60 particles. Usually, however, one is only 
interested in the result on scales of several nanoseconds as 
defined by the interesting frequency range of tens to hun- 
dreds of MHz. On the other hand, it would be useful for 
diagnostic and verification purposes to have a means to 
record the events in the very high time-resolutions of the 
individual pulses. We have therefore implemented both 
possibilities into the simulation program. The different de- 
mands are met by the use of two very different gridding 
strategies. 



Fig. 3. The economic gridding mechanism: when new 
contributions are registered onto existing contributions, 
points are inserted and interpolated only as needed. 



The more efficient and thus favoured strategy for low 
time resolutions is the use of a "simple grid" . It is based on 
an equidistant grid of user-defined resolution (typically ~ 
1 ns). The contributions from a particle's time-series data 
are then binned and thus time-averaged onto the grid. 
The gridding is dynamic in the sense that data points are 
inserted automatically when needed. This saves memory 
in comparison with an ordinary equidistant grid. Still, this 
strategy does not scale well to high time-resolutions. 

To counter the high memory demands for high time 
resolutions, we have implemented an "economic grid" 
which describes the pulse-shapes using only a minimum 
number of points. It is based on an underlying equidistant 
time-grid which limits the time-resolution to a pre-defined 
maximum value. The time-series of an individual particle 
is interpolated to the underlying grid positions and then 
incorporated into the pre-existing time-series, correctly in- 
terpolating any contributions that were already present. 
Points are, however, only inserted in the grid where nec- 
essary. See Fig. 01 for an illustration of the algorithm. 

As another major advantage, the availability of these 
two very different gridding strategies allows an indepen- 
dent cross-check of their implementation. 

4.4. Sequentialised and parallelised calculation 

The program is designed such that calculations can easily 
be sequentialised or parallelised to cut down on memory 
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requirements or take advantage of multiple computers, re- 
spectively. 

Memory usage increases with the number of ground- 
bins to be calculated. To facilitate the use of low-memory 
machines, the calculation of ground-bins can be sequen- 
tialised: the emission from the complete shower is calcu- 
lated for only a subgroup of ground-bins, and only after 
the calculation has finished and the results have been writ- 
ten to disk, the next group of ground-bins is evaluated. 
This efficiently decreases the memory-usage as compared 
to concurrent evaluation of all ground-bins. The overhead 
introduced due to the necessary multiple creation of parti- 
cle distributions is negligible for most combinations of pa- 
rameters. Manually specifying the random seed value for 
the random number generator guarantees identical par- 
ticle distributions in the different calculation segments. 
Allowing for different seed values in the calculations, on 
the other hand, provides a consistency check based on the 
underlying symmetries in the emission pattern as the dif- 
ferent ground-bins are calculated based on independent 
sets of random numbers. 

Similarly, different subgroups of the desired ground- 
area can be calculated in parallel on different computers, 
yielding up to a linear decrease of net computation time. 

4.5. Automatic ground-bin inactivation 

As the simple discrete cutting off of regions outside the 
~ 7 _1 ground-trace region described in section intro- 
duces errors that are too big if one is interested in regions 
of a few hundred metres distance to the shower core, we 
developed a more sophisticated means of cutting down on 
computation time. 

Since most of the particles are distributed in the inner- 
most centre-region of the shower, the ground-bins close to 
the centre-region receive a high number of strong contribu- 
tions, whereas the far-away regions only receive a smaller 
number of not-so-strong contributions. While the centre- 
regions might already have reached sufficient precision af- 
ter a certain number of particles, the outer regions might 
still not have reached the desired statistical precision, af- 
fording a calculation with an even higher number of par- 
ticles. 

The computing time can be distributed in a much more 
efficient way by an on-the-fly inactivation of ground-bins 
that have reached the desired precision. To accomplish 
this, the program evaluates the shower emission in steps 
of (typically) 10,000 particles at a time. After each of these 
steps it compares the (smoothed) time-series derived for 
a specific ground-bin up to that stage with the results for 
that ground-bin at the previous step. Once the relative 
changes fall under a pre-defined limit for a user-defined 
number of steps in a row, the corresponding ground-bin 
is marked as "inactive" and is not evaluated any further. 
The computing time is thus effectively redistributed to the 
outer bins. 



200 




20 40 60 80 100 120 

t [ns] 



Fig. 4. Time-dependence of the raw pulses originating 
from the shower maximum as observed by an observer at 
(from left to right) 20 m, 140 m and 460 m to the north 
from the shower centre. 

The fact that the inactivation sequence should prop- 
agate from the inner to the outer regions is used as a 
consistency check for the procedure. Another advantage 
is that the user can directly specify a desired precision 
rather than having to estimate the number of particles 
needed to reach adequate results. 

5. Data output and reduction 

To facilitate the understanding of the following discus- 
sions, we give a short overview of the raw data that the 
simulations produce and the kind of data "reduction" that 
we apply to visualise the results. 

5.1. Raw data 

The Monte Carlo code tracks the individual particles 
and calculates the associated (vectorial) electromagnetic 
pulses a specific observer (i.e., a specific ground-bin) re- 
ceives. The individual pulses are extremely short, of or- 
der a few 10 -11 s (cf. Fig. [21 . Superposition of all the 
individual pulses yields the raw output of the program: 
one data file per ground-bin stating the time-dependence 
of the north-south, east-west and vertical polarisation 
components of the electric field. Fig. 0] shows the total 
field strength of the raw pulses at different distances to 
the north from the centre of a shower slice consisting of 
10 s particles in the shower maximum. Due to the parti- 
cle's spread in position and time, the pulses are consider- 
ably broader than the individual pulses, of order tens of 
nanoseconds. A helper application is then used to process 
these raw data and reduce them to the desired physical 
quantities. 

5.2. Spectral filtering 

In a first step, the time-series data of the electric field vec- 
tor are Fourier-transformed, yielding the associated spec- 
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Fig. 5. Spectra of pulses originating from the shower 
maximum for observers at (from top to bottom) 20 m, 
140 m, 340 m and 740 m distance to the north of the 
shower centre. 
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Fig. 6. Comparison of the east-west component of a raw 
pulse (solid), a pulse smoothed with a 40-160 MHz ide- 
alised rectangle filter (long dashed) and a 42.5-77.5 MHz 
filter as used in LOPES (short dashed) for emission from 
the shower maximum. 



tra depicted in Fig. |SJ Due to coherence losses caused 
by interference effects, the spectra fall off steeply towards 
high frequencies. At a certain frequency, dependent on the 
distance from the shower centre, the field strength reaches 
a first interference minimum followed by a rapid series of 
alternating maxima and minima in the incoherent regime. 
Insufficient sampling of these extrema yields the unphysi- 
cally seeming features seen here at high frequencies. In the 
emission from a real air shower such features are unlikely 
to exist as the inhomogeneities present in an air shower, 
but not taken into account in the analytic parametri- 
sations, destroy the pronounced extrema. Calculation of 
the emission in this region therefore requires a more de- 
tailed air shower model, e.g. by interfacing of our code to 
CORSIKA. 

Any concrete experiment will have a finite frequency 
bandwidth. Thus, we filter the spectra to infer the ac- 
tual pulses that the experiment will register. As is obvious 
from the spectra, most of the power resides at the low fre- 
quencies. For this reason the amplitude drops significantly 
when filtering frequencies below, e.g., 40 MHz as seen in 
Fig.0 

Here we have used the same idealised 40-160 MHz 
filte r that we app l ied in the theoretical calculations 
of iHueee fc Falckei (|2003^l as well as the actual 42.5- 
77.5 MHz filter that is used in LOPES. The acausality 
of the pulse filtered with the idealised rectangle filter il- 
lustrates that such a filter is an unphysical concept. A 
physical filter does not show acausal behaviour despite its 
steep edges as it delays the signal appropriately. This is 
well visible for the LOPES filter. 



pulse maximum amplitude" can then be visualised in a 
number of ways, e.g., as surface plots, contour plots or 
cuts in specific directions (see the following sections for 
examples) . The absolute time associated to the maximum 
pulse amplitude additionally yields information about the 
curvature of the radio wave front. 

One subtlety involved with this procedure is the noise 
levels present at high radial distances. As can be seen in 
Fig. the pulses get significantly broader as one goes 
to higher distances. (This effect gets even stronger for 
fully integrated showers.) At distances of several hun- 
dred metres, the pulses are so broad that a filter clip- 
ping frequencies below ~40 MHz actually resolves the 
pulses out. Consequently, the pulse amplitude drops to 
very low values comparable to those introduced by the 
higher-frequency numerical noise associated with the very 
short individual particle pulses. Taking the maximum am- 
plitude as a measure for the emission strength then might 
no longer constitute a useful procedure. Calculating the 
time-integral of the field strength or the received power 
would be a better approach in these cases. 

For the diagnostics of the employed algorithms in sec- 
tion El and the analysis of track length effects in section 
17.11 we therefore adopt 800 m as a cutoff distance. For 
the calculation of emission strengths from a shower slice 
and an integrated shower, we limit the plots to distances 
of ^550 m and ~400 m, respectively, to insure that the 
filtered pulse amplitudes give adequate results which are 
directly related to the experimentally relevant signal-to- 
noise levels. 



5.3. Further data processing 

To analyse the radial dependence of the emission strength, 
wc then determine the maximum amplitude of the filtered 
pulse in each ground-bin. The dependence of this "filtered 



6. Consistency checks 

We have studied the output of our MC code very carefully 
to make sure that the calculations are correct. In particu- 
lar, we have made the following consistency checks: 
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Fig. 7. Total field strength of a pulse originating from 
a point source of 10 8 7 = 60 particles in 4 km height, 
observed in the shower centre. Solid: analytic calculation, 
points: calculated with MC code. 
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Fig. 8. East-west polarisation component of a point- 
source with 7 — 60 particles at 4 km height at increasing 
distance to the north from the shower axis. From left to 
right: 5 m, 305 m, 505 m, 705 m, 855 m and 995 m. 



6.1. Individual particle pulses 

Figure [3 shows the comparison between an analytical cal- 
culation and our MC code for a pulse created by a point- 
source consisting of 10 8 particles with 7 = 60 at 4 km 
height, comparable to a pulse that would originate from 
the maximum of a vertical 10 17 eV air shower concentrated 
into a point, as seen from an observer in the shower cen- 
tre. The magnetic field is adopted as horizontal with a 
strength of 0.3 Gauss. The good agreement between the 
MC and analytic results demonstrates that the calculation 
of the particle trajectories and emission contributions is 
implemented correctly. 

Figures [S] and [5] demonstrate how the pulses from the 
same point-source change when the observer moves out- 
wards from the shower centre to the north and east, re- 
spectively. While the pulse amplitude drops quickly when 
one goes to the north, it stays fairly constant as one goes 
to the east. This ef fect was already discus sed under the 
term "reduced 9" in lHuege fc Falckd l)2003|) : The particle 
trajectory bends towards an observer in the cast or west, 
so that he or she still sees the particle with a very small 
angle to the line of sight — just during a different part 
of the trajectory. For the specific geometry chosen in this 
example, i.e. a particle pair starting off vertically down- 
wards in the shower centre, the pulses even get broader 
as one goes outwards to the east, since one only sees an 
(asymmetric) half pulse in the centre and can sec the full 
(symmetric) pulse only at considerable distance. All in all 
the behaviour is exactly as expected. 

Another important consistency check is the depen- 
dence of the individual pulses on the magnetic field 
strength. Figure shows the changes arising when the 
magnetic field is enhanced from 0.3 Gauss to 0.5 Gauss 
(both with 0° inclination, i.e. horizontal). The pulse gets 
stronger, yet at the same time shorter so that the inte- 
gral over E dt stays constant (while the power an observer 
receives from an individual particle, i.e. the integral over 
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Fig. 9. Same as figure[S]at increasing distance to the east 
from the shower axis. 

E 2 dt, thus scales linearly with B), exactly as inferred from 
analytic calculations. This directly leads to an important 
result: The overall radio emission from an air shower can- 
not depend strongly on the magnetic field strength, as it 
is the sum of a great number of individual pulses, the in- 
tegral of each of which is independent of the value of B. 
Likewise, asymmetries between north and south that are 
introduced by a realistically inclined magnetic field cannot 
be very strong. 

The fact that the 0.5 Gauss pulse arrives earlier than 
the 0.3 Gauss pulse for an observer 995 m to the east of 
the centre is explained by the stronger curvature of the 
particle trajectory. 

6.2. Symmetry N-S and E-W 

For a vertical air shower, the emitted radiation pattern 
must have a number of inherent symmetries. In particular, 
the east and west directions are completely equivalent (as 
long as particle track lengths for electrons and positrons 
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Fig. 10. Individual particle pulses from a point-source 
shower as measured by an observer situated 995 m east 
of the shower centre. Solid: B = 0.3 Gauss, dashed: B = 
0.5 Gauss. The time-integral over the pulses is constant. 



are adopted as identical), so that the emission pattern 
must be symmetric in east and west. For a horizontal mag- 
netic field, north and south must also be equal. The fact 
that these intrinsic symmetries arc fulfilled as expected 
provides an important consistency check. Furthermore, 
one can use these symmetries to save computation time 
by calculating only a half- or even quarter-plane on the 
ground and then mirroring the results accordingly. 

6.3. Gridding algorithms 

As explained in section T4. 31 we have implemented a "sim- 
ple grid" for low time-resolution calculations as well as an 
"economic grid" for cases in which the user is interested in 
the full time-resolution associated to the individual parti- 
cle pulses. 

The independence of the two algorithms allows a cross- 
check of the routines. Fig. II II shows the raw pulse as cal- 
culated with the different gridding strategies at different 
time resolutions. The result is consistent between all three 
cases, which demonstrates that both algorithms work well. 
The "simple grid" proves to be very efficient at low time- 
resolutions (typically ~ 1 ns). It is especially robust in 
the sense that it remains stable regardless of the specific 
resolution used. The "economic grid" on the other hand 
has to be set to a time-resolution high enough to resolve 
the individual particle pulses. Especially when particle en- 
ergy distributions arc switched on, introducing very high- 
energy particles with up to 7 = 1000, this strategy quickly 
becomes inefficient. 

Additionally, low time-resolution "simple grid" data 
yields smoother pulses due to the time-averaging over the 
individual particle pulses, which allows precise calcula- 
tions with fewer particles. 
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Fig. 11. Raw pulse from a shower slice calculated with 
different gridding strategies and resolutions. Solid with 
points: simple grid 10~ 9 s, solid black: simple grid 10 -10 s, 
light coloured: economic grid 10 _ 12 s. 
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Fig. 12. Changes introduced by the smart sampling al- 
gorithm in the radial emission pattern of the rectangle- 
filtered maximum pulse amplitude for emission from the 
shower maximum. Thin lines: dense equidistant sampling, 
thick lines: smart sampling; solid: to the north, dashed: to 
the west. 



6.4. Smart trajectory-sampling 

Figure IT21 demonstrates that the errors introduced by the 
smart sampling algorithm are only slight. At the same 
time, the algorithm allows a huge cut-down on compu- 
tation time. The algorithm can, however, optionally be 
switched off for a more precise calculation. 

6.5. Cutting off 7 _1 -cones 

As mentioned earlier, the discrete cutting off of radiation 
contributions outside the ground-trace of a few 7 _1 -cones 
width can decrease the computation time enormously. 
FigurcEH however, demonstrates that this strategy is not 
suitable if one needs precision at the percent-level to be 
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Fig. 13. Changes introduced by the cutting off of regions 
outside a few 7 _1 -cones in the radial emission pattern of 
the frequency-filtered maximum pulse amplitude for emis- 
sion from the shower maximum. Distance from the shower 
centre is to the north. Solid: no cutting, long dashed: cut- 
ting after 8 7 -1 , short dashed: cutting after 4 -f^ 1 
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Fig. 14. Changes introduced to the radial emission pat- 
tern of the frequency-filtered maximum pulse amplitude 
for emission from the shower maximum by the automatic 
bin inactivation algorithm. Thin lines: no automatic bin 
inactivation, thick lines: automatic bin inactivation; solid: 
to the north, dashed: to the west. 



able to describe the emission pattern out to distances of 
several hundred metres from the shower centre. 

The discrete cutting introduces "breaks" in the radial 
emission pattern exactly at the positions corresponding 
to the cutoff, i.e. at ks 530 m in case of 8 7~ 1 -cones and 
« 270 m in case of 4 7 _1 -cones for 7 = 60 particles at 4 km 
height. The effect is less strong in the east- west direction, 
but overall this strategy is disqualified for high-precision 
calculations. 

The problems are resolved by the more sophisticated 
on-thc-fly inactivation of ground-bins. 

6.6. Automatic ground-bin inactivation 

Fig. [H] demonstrates the stability of the automatic 
ground-bin inactivation algorithm. The calculations are 
as precise as those with a fixed high number of particles 
for all ground-bins, yet at the same time allow a much 
more efficient use of computing time and the definition of 
a user-specified precision goal. The deactivation sequence 
propagates from the inside to the outside as expected (see 
Fig. ll5J) . It is also no surprise that the bins in the far east 
and west require the most computing time, as these are the 
regions that are influenced most by edge effects from tra- 
jectory cutoffs, different trajectory curvatures and the like. 
The automatic ground-bin inactivation strategy therefore 
turns out to be a very powerful and self-consistent tech- 
nique to conduct extensive simulations with high preci- 
sion. 

7. Emission from a shower slice 

Similar to the ana lytical calculations described in 
iHuege fc Falckel l|2003h . we first take a look at the emis- 
sion from a single "slice" of the air shower. Throughout 
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Fig. 15. Automatic ground-bin inactivation sequence. 
Darker bins are set inactive later than lighter bins. The se- 
quence propagates from the inside to the outside. The pat- 
tern is east-west and north-south symmetric as expected 
for a vertical shower and a horizontal magnetic field. 



this section, we consider only the maximum of a vertical 
air shower induced by a 10 17 eV primary particle, consist- 
ing of 10 s charged particles at a height of 4 km. In a step 
by step analysis, we increase the complexity of the particle 
distributions and evaluate the changes introduced in the 
simulation results. 
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Fig. 16. Radial dependence of the maximum rectangle- 
filtered pulse amplitude for emission from the shower 
maximum in the north (solid) and west (dashed) direc- 
tion in case of constant and long particle trajectories 
(A = 100 g cm" 2 , no edge effects at high distances from 
the shower centre). 



7.1. Trajectory length effects 

In a first step, we adopt the geomagnetic field parallel to 
the ground with a strength of 0.3 Gauss, as it is present at 
the equator. This is the same configuration that we used in 
the theoretical calculations. For the moment, we consider 
the simplified case of monoenergetic 7 = 60 particles. 

The theore t ical calculations carried out in 
iHuege k Falckel (2003) were based on an analytical 
derivation of the spectra of individual particles on 
circular trajectories. This derivation makes the implicit 
assumption that the trajectory is always symmetric with 
respect to the point in which the minimum angle 9 to the 
observer's line of sight is reached. In other words, edge 
effects arising from the cutting off of the finite trajectories 
are not taken into account. 

These edge effects, however, turn out to significantly 
shape the radial emission pattern of the radiation. Let 
us first consider the case of trajectories which are long 
enough to not produce significant edge effects for observers 
far away from the shower centre by adopting a trajec- 
tory length of 100 g cm -2 . The result is depicted in Fig. 
1161 As expected, it is very similar to the theoretical pre- 
diction for the " reduce d 9" case presented in Fig. 14 of 
iHueee k Falcke] l|2003|) . In this scenario, one would ex- 
pect a significant asymmetry between the north-south and 
east-west direction, which would obviously be a very use- 
ful observable. 

What happens, however, if one adopts more realistic 
particle track lengths? Let us first consider a constant tra- 
jectory length of 40 g cm" 2 , which is approximately the 
free path length (equal to one radiatio n length) of elec- 
trons and positrons in air llAllanl . ll97l|) . In this scenario, 
edge effects strongly shape the emission at high distances 
in the east- west direction as shown in Fig. El An unusual 
"kink" appears at ^540 m. This is not a numerical glitch, 
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Fig. 17. Same as Fig.Elfor A = 40 g cm -2 . See text for 
explanation of the "kink" at ~540 m. 
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Fig. 18. Time-dependence of the raw pulses originating 
from the shower maximum as observed by an observer 
at (from left to right) 460 m, 500 m, 540 m, 580 m and 
620 m to the west from the shower centre. See text for 
explanation of the "polarity change" . 

but an interference effect arising when the observer stops 
to see the main (positive) peak of the individual particle 
pulses due to the edge effects. He then only receives the 
initial (negative) contribution of the electric field pulse. 
This effectively causes a polarity change in the raw pulses 
as shown in Fig. 1181 accompanied by a temporary drop in 
the filtered pulse amplitude. 

Finally, we change to a realistic exponential distribu- 
tion of track lengths with a mean of A = 40 g cm" 2 , 



p(X) = p exp 



X 

T 



(5) 



As can be seen in Fig. 1191 the asymmetry between north- 
south and east-west direction is now washed out up to high 
distances. (In fact, it will be washed out almost completely 
once the integration over the shower evolution is taken into 
account.) 

Apart from the regions far from the shower centre, 
edge effects also occur in the centre region due to the 
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Fig. 19. North (solid) vs. west (dashed) asymmetry in 
the maximum filtered pulse amplitude for emission from 
the shower maximum in case of statistically distributed 
track- lengths. The asymmetry is washed out up to high 
distances. 

instantaneous starting of the trajectories. These effect are 
discussed in section I7~H 

The significance of the trajectory length effects already 
illustrates the importance of adopting as realistic proper- 
ties for the particle distributions as possible, a goal that 
could not be reached with analytic calculations alone. We 
retain the realistic statistical distribution of track lengths 
for the following discussions of the magnetic field and en- 
ergy distribution effects. 

7.2. Magnetic field dependence 

We now make the change from an equatorial 0.3 Gauss 
horizontal magnetic field to the 0.5 Gauss 70° inclined 
magnetic field present in central Europe. As explained in 
Sec. 16.11 the influence of the magnetic field on the inte- 
grated shower pulse should not be very significant. This 
is confirmed by Fig- EDI Although the projected magnetic 
field strength drops from 0.3 Gauss to 0.17 Gauss, the 
emission pattern only changes very slightly. In particular, 
the amplitude level stays almost constant. Interestingly, 
the north-south vs. east-west asymmetry is washed out 
even further. 

On the other hand, the inclination of the magnetic field 
breaks the intrinsic north-south symmetry of the emission 
pattern. This is visible in the north-south and east-west 
polarisation components of the contour plots shown in Fig. 
1211 Similarly, the symmetry of the ground-bin inactiva- 
tion sequence is broken as shown in Fig. |23 The effect 
is, however, only weak and there is still no asymmetry if 
one considers the total electric field strength rather than 
a specific polarisation direction. 

Note that there is no significant emission in the centre 
region of the north-south polarisation component (second 
column of Fig. I21JI . This demonstrates that the numer- 
ical cancellation of north-south polarised radiation com- 
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Fig. 20. Changes to the north (solid) and west (dashed) 
radial emission patterns for emission from the shower max- 
imum when going from a 0.3 Gauss horizontal magnetic 
field (thin lines) to a 70° inclined 0.5 Gauss magnetic field 
(thick lines). 
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Fig. 22. Automatic ground-bin inactivation sequence in 
case of 70° inclined 0.5 Gauss magnetic field. The north- 
south symmetry is broken as expected, cf. Fig. 1151 

ponents from the elect rons and positrons (labelled A±_ in 
iHueee k Falckel i|2003F) ') indeed works correctly in the MC 
code and that it was justified to neglect this component 
in the earlier theoretical calculations at least for the pri- 
marily important centre region. 

In summary, the magnetic field effects are weak, but 
non-trivial. 

7.3. Energy distribution effects 

In the theoretical calculations, a change from monoen- 
ergetic 7 = 60 electrons to a broken power-law peak- 
ing at 7 = 60 introduced only minor changes, namely a 
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Fig. 21. Contour plots of the 40-160 MHz rectangle-filtered maximum pulse amplitude for emission from the shower 
maximum in case of a horizontal 0.3 Gauss magnetic field (upper panel) and a 70° inclined 0.5 Gauss magnetic field 
(lower panel) and 7 = 60 particles. Contour levels are 5 /iV m _1 apart. From left to right: total electric field strength, 
north-south polarisation component, east- west polarisation component, vertical polarisation component. 



slight redistribution of radiation from medium distances 
to the centre region (due to the high-energy particles) 
and high distances (due to the low-energy particles). The 
changes introduced by the energy distribution are bound 
to be more complex in the MC simulations, as they are 
heavily intertwined with the edge effects discussed before. 
Nonetheless, a similar redistribution of radiation from the 
medium scales to the innermost centre region is observable 
in the MC simulations as shown in Fig. 1231 

The drop visible in the west direction at ^200 m arises 
due to a transition from a dominating east-west polari- 
sation component to a dominating north-south polarisa- 
tion component in combination with some resolving out of 
pulses due to the filter bandwidth used. (Fig. [2] demon- 
strates that the north-south and east-west polarisation 
components become comparable at these distances.) 

It should be pointed out, however, that the strength 
of the effects introduced by the choice of a specific energy 
distribution seem misleadingly strong when one consid- 
ers the emission from a single shower slice alone. Once 
the integration over the shower evolution as a whole is 
performed, most effects (including the drop at ^200 m) 
are again washed out almost completely and the influence 
of the specific choice of energy distribution becomes very 
weak. For the moment, we therefore continue to use the 
broken power-law distribution that we adopted in the the- 
oretical calculations rather than implement ing a more re- 
alistic distribution such as the one given bv lNerling et alJ 
ll2003h . 

The result that we have reached so far is the emis- 
sion from the maximum of a 10 17 eV air shower consisting 
of 10 s electrons and positrons at a height of 4 km, tak- 
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Fig. 23. Changes introduced when switching from mo- 
noenergetic 7 = 60 particles (thin lines) to a broken 
power-law peaking at 7 = 60 (thick lines) for emission 
from the shower maximum. Solid: to the north, dashed: to 
the west. See text for explanation of the drop at ~200 m. 

ing into account adequate spatial, energy and trajectory 
length distributions and a magnetic field as present in cen- 
tral Europe. It is illustrated once more as a surface plot 
in Fig. El 

7.4. Comparison with theoretical calculations 

As shown in Fig. the individual particle pulses calcu- 
lated with our MC code are consistent with the results of 
the analytic calculations. The analytic pulses are symmet- 
ric because the particle trajectories are implicitly adopted 
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Fig. 24. Pattern of the maximum filtered electric field 
amplitude (in /iV m _1 ) emitted by the maximum of a 
10 17 eV air shower consisting of 10 s particles at a height 
of 4 km for a 0.5 Gauss 70° inclined magnetic field, a 
broken power-law particle energy distribution and a sta- 
tistical distribution of track lengths. 

symmetric to the point of smallest angle to the line of 
sight — the analytic calculations do not consider edge ef- 
fects associated with the finite lengths of the trajectories. 
When one takes into account the cutting off of the tra- 
jectories correctly in the MC calculations, an observer in 
the centre region only sees half of the symmetric pulse for 
each individual particle (cf. Fig. El> • 

In section T6. II we discussed that the time integral over 
the individual particle pulses is the quantity relevant for 
the overall pulse as integrated over the shower as a whole. 
For simple geometries such as a point source or a line 
charge, the change from the symmetric to the half pulses 
should therefore produce a drop in the overall pulse am- 
plitude (and the overall pulse spectrum) by a factor ~ 2. 
(For more complex shower geometries or at higher dis- 
tances from the shower centre the changes introduced by 
the edge effects are non-trivial but still important as al- 
ready discussed in section mi - ) 

Fig. 1251 shows the direct comparison between the ana- 
lytic and MC simulated spectra emitted by a slice of mo- 
noenergetic 7 = 60 particles for the case of a horizontal 
0.3 Gauss magnetic field and long (constant 100 g cm -2 ), 
but not symmetric, particle trajectories. (Because the 
north-south polarisation component was neglected in the 
analytic calculations, we hereafter directly compare the 
east-west polarisation component.) This scenario allows 
a very direct comparison between the analytic and MC 
calculations, the only major difference being the edge ef- 
fects introduced due to the non-symmetric trajectories. 
Consequently, the analytic spectra lie a factor ~ 2 above 
the MC results. Scaling down the analytic results by a fac- 
tor of two shows indeed a very good agreement between 
the analytic and MC results as seen in Fig. The radial 
dependence of the emission pattern, shown in Fig. 1271 also 
shows good agreement between the MC and analytic cal- 
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Fig. 25. Spectra emitted by the maximum of a 10 17 eV 
vertical air shower consisting of 10 s particles with 7 = 
60 at a height of 4 km, long (constant 100 g cm -2 ), 
but not symmetric, particle trajectories and horizontal 
0. 3 Gauss magnetic field. Thin lines: analytic calculations 
of iHuege fc Falckel l)2003j) . thick lines: these MC simula- 
tions. Solid: 20 m, dashed: 100 m, dotted: 260 m to north 
from shower centre. 
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Fig. 26. Same as Fig. ESI but scaling down the analytic 
results by a factor of two. 

culations when one accounts for the systematic factor of 
two. 

In a next step, we switch on the statistical distribution 
of trajectory lengths with a realistic mean free path length 
of 40 g cm -2 , adopt again a broken power-law distribution 
of particle energies and change to the realistic 70° inclined 
0.5 Gauss magnetic field present in Central Europe. In 
the analytic calculations, we switch on the broken power- 
law distribution, but cannot take into account the track 
length effects or the inclined magnetic field. Obviously, the 
results produced by the MC simulations in this scenario 
therefore cannot be expected to reproduce the analytic 
results equally well. 

A direct comparison between the spectral dependences 
predicted by the analytic and MC calculations is shown 
in Fig. EHl keeping the down-scaling of the analytic re- 
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Fig. 27. Radial dependence of the emission for the same 
scenario as in F i g. 1251 Thin lines: analytic calculations of 
iHueee fc Falcke] (2003) scaled down by a factor of two, 
thick lines: these MC simulations. Solid: v = 50 MHz, 
dashed: v = 75 MHz, dotted: v = 100 MHz. Distance is 
to north from shower centre. 
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Fig. 28. Spectra emitted by the maximum of a 10 17 eV 
vertical air shower consisting of 10 8 particles with broken 
power-law energy distribution at a height of 4 km, sta- 
tistically distributed particle trajectories with 40 g cm -2 
mean path length and 70° inclined 0.5 Gauss magnetic 
field. Thin lines: analytic calculations of iHueee fe FalckH 
( 2003) scaled down by a factor of two, thick lines: these 
MC simulations. Solid: 20 m, dashed: 100 m, dotted: 260 m 
to north from shower centre. 

suits by a factor of two to compensate for the symmetric 
trajectories. As discussed earlier, switching on the energy 
distribution redistributes flux from medium scales to the 
centre region. Correspondingly, the spectrum at 20 m dis- 
tance shifts to higher amplitudes. The effect is stronger 
in the MC calculations than in the analytics. The MC 
100 m spectrum fits well with the analytic results, whereas 
the MC 260 m spectrum lies at lower amplitudes than in 
the analytics: The radial dependence now falls off much 
steeper in the MC as compared to the analytic results, as 
is also visible in Fig. |2H1 Overall, the radial dependence 
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Fig. 29. Radial dependence of the emission for the same 
scenario as in F i g. I28I Thin lines: analytic calculations of 
IHueee fc Falcke] (2003) scaled down by a factor of two, 
thick lines: these MC simulations. Solid: v — 50 MHz, 
dashed: v = 75 MHz, dotted: v = 100 MHz. Distance is 
to north from shower centre. 

follows a very much exponential decay. (We will carry out 
a detailed analysis of the functional form of different prop- 
erties of the radio emission in a later paper.) It does not 
exhibit the prominent plateaus visible in the centre re- 
gions of the analytic calculations. Considering the much 
higher precision of the MC simulations and the neglect of 
edge effects and statistical trajectory lengths in the an- 
alytic calculations, however, we consider the agreement 
quite acceptable. 

8. Emission from an integrated shower 

After having analysed the emission from a shower slice, the 
next step in our analysis now is to perform the integration 
over the air shower as a whole. 

8.1. Integration over shower evolution 

In the theoretical calculations we performed in 
iHuege fc Falcke] fj2003), the integration over the shower 
evolution was carried out in a somewhat simplified 
way: The shower evolution was discretised into slices 
of independent generations of particles, spaced apart 
by one radiation length each. The overall emission was 
then superposed as the sum of the radiation from all 
these slices. Although this strategy should allow a good 
estimate of the emission from the complete shower, it has 
at least two problems: First, the total number of particles 
in the shower — and thus the total field amplitude — is 
directly influenced by the scale introduced through the 
spacing of the slices. A denser or wider spacing directly 
leads to higher or lower emission levels, respectively. 
Although the radiation length is the logical choice for this 
scale, a "scale-free" approach would be a better choice. 
Second, our theoretical calculations, strictly speaking, 
are only valid in the far field. Consequently, the emission 
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from slices close to the ground, especially for high-energy 
showers, cannot be taken into account with the desired 
precision. 

Both these pitfalls no longer pose a problem in the 
MC simulations: No far- field approximations at all were 
applied in the MC calculations, and the continuous evo- 
lution of the shower is correctly taken into account in the 
creation of particles according to the corresponding prob- 
ability distribution function. 

For the show e r prof ile we use the Greisen parametri- 
sation llGreisenl Il960l) that we already adopted in 
iHuege fc FaickelfcOO.Sl) : 



N(s) 



0.31 exp 



X m 2-3 Ins 



X 3/s-l 



V X m / Xq 



(6) 



where the (theoretical) position of the shower maximum 
X m is given by 



X m — X In (E p / E CI it) 



(J) 



Xq = 36.7 g cm -2 denotes the electron radiation length in 
air, iS cr it = 86 MeV corresponds to the threshold energy 
where ionisation losses equal radiation losses for electrons 
moving in air, and E p specifies the primary particle energy. 
Equation (jTJ predicts the position of the shower maximum 
for a purely electromagnetic cascade, in which the shower 
age as a function of atmospheric depth then corresponds 
to 



s(X) 



3X 



X + 2X n 



(8) 



Obviously, this parametrisation of the shower age has to 
be modified to adequately describe the evolution of the 
hadronic air showers in our MC simulations. We choose 
to manually set the depth of the shower maximum to an 
empirical value A" m . c as a function of primary particle en- 
ergy and shower inclination, e.g. X m . c ~ 630 g cm -2 for 
our typical 10 17 eV vertical air shower. The shower age is 
then adopted as 



s(X) 



3X 



X + 2X„ 



(9) 



whereas eq. © is left unchanged, i.e. retaining the theo- 
retically motivated value for X m . 2 

Eq. ijHJ denotes the integrated number of electronic 
particles that a detector positioned at atmospheric depth 
X measures as the shower sweeps through it. This number 
is not equal to the number of particles I(X) that are "in- 
jected" at that atmospheric depth, the quantity we need 
to describe the probability distribution function for the 
creation of particles. The two quantities are directly re- 
lated via the path length distribution of the particles. For 

2 A more realistic set of parameters for application of the 
Greisen function to hadronic showers in th e energy range be- 
tween 10 17 and 10 18 eV was established in lAbu-Zavvad et alJ 
(2001). For the moment, however, we retain the parametrisa- 
tion as stated above to allow a better comparison with our 
earlier theoretical calculations. 
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Fig. 30. Comparison of I(X) (solid) and N(X) (dashed) 
as a function of atmospheric depth X for a vertical 10 17 eV 
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Fig. 31. Trace of the trajectories in a complete 10 17 eV 
air shower. 

an exponential path length distribution with mean free 
path length A as given in eq. J5J, the injection function is 
given by 



I(X) 



dN 
dX 



N 

y 



(10) 



As demonstrated in Fig.[3U| I(X) closely follows the form 
of N(X) with an offset of ~ A to lower X values, i.e. 



I(X) &N(X + \)/\. 



(11) 



The shower evolution is thus taken into account in a 
continuous and consistent way by random creation of par- 
ticles with probabilities according to eq. i|10|) . Fig. [^illus- 
trates the shower evolution through the particle trajecto- 
ries that are followed during the simulation of our 10 17 eV 
vertical air shower. 

8.2. Integrated shower results 

As in the theoretical calculations, the integration over the 
shower evolution has two main effects, visible in Fig. 1321 
First, the emission level is boosted significantly. This di- 
rectly shows that the emission is not described sufficiently 
by just taking into account the shower maximum. The 
second major effect is a steepening of the radial emission 
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Fig. 32. Effects introduced in the radial dependence of 
the 40-160 MHz rectangle-filtered pulse amplitude by the 
integration over the shower evolution. Thick lines: inte- 
grated 10 17 eV vertical shower with broken power-law 
particle energy distribution, statistical track length dis- 
tribution with A = 40 g cm -2 and 70° inclined 0.5 Gauss 
magnetic field, thin lines: only shower maximum in same 
scenario. Solid: to the north, dashed: to the west. 

pattern due to the amplification of coherence losses. It is 
mainly the centre region which receives significant addi- 
tional radiation. The steepness of the radial dependence is 
also illustrated by the strong drop in the maximum ampli- 
tude of the filtered pulses when one goes to even moderate 
distances of 260 m as shown in Fig. 02I 3 

Another important effect is the further fading away of 
sharp features and the asymmetries associated with the 
geomagnetic field in the emission pattern, as illustrated 
by Fig. |21 While some of these features were still quite 
prominent in the emission from a single shower slice, they 
more or less vanish completely as soon as the integration 
is performed. Similarly, the effects of individual changes to 
the particle distributions, such as the specific choice of an 
energy distribution, are therefore much weaker than they 
are in case of a single shower slice. 

The polarisation characteristics of the integrated 
shower are very similar to those of an individual slice. 
Again, the emission is highly polarised in the east- west 
direction, except for some regions at medium distances 
of ^200-400 m as illustrated in Fig. 021 and, in a more 
quantitative way, Fig. 1351 

Figure OH shows a direct comparison of the radial 
profiles of the particle density (as given by the NKG- 
parametrisation) and the filtered radio pulse amplitude. 
The particle density falls off much more steeply than the 
radio signal in the central ~200-250 metres. Further out, 
the slope becomes comparable for radio amplitude and 
particle density. 

Th e integrated pulse shown in Fig. 18 of lHuee c & Falcke 
(2003) was by mistake calculated for a 40-160 MHz rectangle 
filter. Using the intended 42.5-77.5 MHz rectangle filter, the 
amplitude drops from ~ 1500 //V m -1 to ~ 1100 m _1 . 
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Fig. 33. Pulses in the east-west polarisation component 
after 40-160 MHz rectangle-filtering for a shower as de- 
scribed in Fig. |23 Solid: in the shower centre, long dashed: 
100 m to north of centre, short dashed: 260 m to north of 
centre. 
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Fig. 35. Radial dependence of E(R, u>) for different polar- 
isation components at v = 55 MHz for the same scenario 
as in Fig. |2U Solid: east-west polarisation to the north 
from centre, long-dashed: north-south polarisation to the 
north from centre, short-dashed: east-west polarisation to 
the north-west from centre, dotted: north-south polarisa- 
tion to the north-west from centre. 

8.3. Comparison with theoretical calculations 

We now compare the results of our MC simulations of a 
fully integrated 10 17 eV vertical air shower with the the - 
oretical calculations performed in iHueere fc Falcke] ll2003l) . 
Fig. shows the spectral dependence of the emission in 
direct comparison. The MC results again produce some- 
what lower levels of radiation. Scaling down the theoretical 
results by the systematic factor of two introduced in sec- 
tion l7.4l the agreement is much better, as shown in Fig. 1381 
Considering the huge differences in the two calculations, 
the a greement is quite remar kable. 

In iHuege fc Fakkl {2003) we compared the theoreti- 
cal results with the available historical data. As discussed 
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Fig. 34. Contour plots of the 40-160 MHz rectangle-filtered pulse amplitude for the same scenario as described in 
Fig. Contour levels are 20 /iV m -1 apart. From left to right: total electric field strength, north-south polarisation 
component, east- west polarisation component, vertical polarisation component. 
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Fig. 36. Comparison of the filtered radio pulse ampli- 
tude (solid) and the particle density as given by the NKG- 
parametrisation (dashed) as a function of radial distance 
from the shower centre. Values were normalised to unity 
at r = 60 m. 

there in detail, the absolute values of the historical exper- 
imental data are very uncertain and largely discrepant be- 
tween the different groups. Additional uncertainty arises 
from ambiguities in the exact definition of the historical 
values denoted as e„ and their conversion to the theoretical 
quantity \E(R, lo)\, which we performed via the relation 

/ 1 no 

e v = J \E(R,u)\ fn6A\E(R,w)\. (12) 

V 7T 

To deal with the discrepa nt sets of da ta we decided to 
take the well documented lAllanl l)l97l|) data as ou r ref- 
erence and res cale the spectral data of IPrahl l)l97l|) and 
ISpencerl l|l969f) to be consistent with the Allan data at 
v = 55 MHz. In this work, we use the identically rescaled 
data for comparison with our new MC spectra. While the 
absolute values of the spectral data are therefore some- 
what arbitrary, they at least allow an evaluation of the 
qualitative spectral dependence. The spectral data are 
over-plotted in Fig. |2H1 The agreement with the spectral 
dependence in the centre region of the shower is very good. 

In Fig. [2HI we compare the radial dependence of the 
emission as calculated by our MC code with the earlier 
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Fig. 37. Spectra emitted by a complete 10 eV vertical 
air shower with maximum at 4 km height, broken power- 
law particle energy distribution, statistically distributed 
particle trajectories with 40 g cm -2 mean path length and 
70° inclined 0. 5 Gauss magnetic field. Thin lines: analytic 
calculations of iHuege fc Falcke] £2003), 

thick lines: these 

MC simulations. Solid: shower centre, dashed: 100 m, dot- 
ted: 250 m to north from shower centre. 

theoretical results and the lAllanl ljl97l[) data. The scaled- 
down theoretical results again show good agreement with 
the MC results in the centre region. Overall, the MC pre- 
dicts a somewhat steeper decrease of the emission levels to 
higher distances. It was, however, expected that the the- 
oretical calculations overestimate the coherence and thus 
emission levels at high distances. The absolute level of the 
emission also agrees with the Allan data within the un- 
certainties. 

9. Discussion 

With the design and implementation of our MC simula- 
tion, we have taken the logical next step in the devel- 
opment of our geosynchrotron radiation model for radio 
emission from cosmic ray air showers. The MC technique 
provides an independent cross-check on our earlier theo- 
retical works and allows us to model the emission from the 
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Fig. 38. Same as Fig. &\ but scalin g down the a nalytic 
res ults by a f actor of two. Data from IPrabl l)l97l[) (gray) 
and lSpencerl lll969t) (b lack) were rescaled to be consistent 
with the lAllanl(|l97l[) data at v = 55 MHz. 
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Fig. 39. Radial dependence of the emission at v = 
55 MHz for the same scenario as in F i g. |37l Dashed: an- 
alytic calculations of Huege & Falcke (2003J). Solid: ana- 
lytic calculations o flHuege fc Falckd l)2003|) scaled down by 
a factor of two (thin line) in comp arison with these MC 
simulations (thick line). Data from lAllanl 1 197lj) . Distance 
is to north from shower centre. 



air shower with much higher precision based on a much 
more realistic air shower model. 

To make these simulations feasible on standard "off- 
the-shelf" computer hardware, we conceived and imple- 
mented a number of intelligent concepts. In particular, 
the smart sampling of the particle trajectories in conjunc- 
tion with the automatic inactivation of ground-bins and 
the adaptive collection of the time-series data provide the 
necessary cut-down on computation time compared to a 
simple "brute-force" approach. At the same time, high 
precision is retained in the results as demonstrated by the 
detailed consistency checks presented in section|H| The cal- 
culation of a typical vertical 10 17 eV air shower with real- 
istic shower properties and 25, 000, 000 particles combined 
with automatic ground-bin inactivation at very high preci- 



sion takes about 400 seconds per ground-bin on a standard 
PC of the 1.6 GHz class. A very detailed calculation with 
800 ground-bins (25 radial, 32 azimuthal), parallelised on 
4 PCs, thus takes about one day. 

A major new result of our MC simulations is the 
profound influence of edge effects arising from the finite 
lengths of the particle trajectories. These could not be 
taken into account in the theoretical calculations. In com- 
bination with the dependence of the individual particle 
emission on the magnetic field strength, the edge effects 
associated with the statistical distribution of track lengths 
wash out the asymmetries originally introduced into the 
total field strength emission pattern by the geomagnetic 
field to a high degree. Once the integration over the shower 
evolution as a whole is carried out, the asymmetry is gone 
completely. The absence of asymmetries in the total field 
strength emission pattern from an integrated shower is 
somewhat unfortunate, as a prominent asymmetry would 
have allowed an easy, yet unambiguous test of the geo- 
magnetic emission mechanism via the statistics of radio 
pulse total field strengths alone. 

The decomposition of the electric field into north- 
south, east-west and vertical polarisation components, 
however, shows that the emission is indeed highly po- 
larised in the direction perpendicular to the shower axis 
and the magnetic field direction, as predicted by our the- 
oretical calculations. Polarisation-dependent radio mea- 
surements such as the ones carried out by LOPES could 
therefore still unambiguously establish the geomagnetic 
emission mechanism. 

The MC results show good consistency with the theo- 
retical spectra and radial dependences — apart from a sys- 
tematic factor of two in emission strength which is plausi- 
ble considering the implicit assumption of symmetric tra- 
jectories in the analytic calculations. Such good agreement 
between the theoretical and MC calculations is remarkable 
considering the inherent differences in the two approaches. 
In particular, the integration over the shower evolution as 
a whole is carried out in a much more sophisticated way in 
the MC simulations as compared to the analytical works. 
As mentioned above, the polarisation characteristics of the 
emission are exactly as inferred in the theoretical calcula- 
tions. 

The spectra and radial dependences predicted by our 
MC co de also agree wel l with the hi storic al data o f lAllanl 
l)l97ll) and the data of lPrahl l|l97lh and ISnencerl |l969h 
scaled to the absolute level of the Allan data. The nec- 
essary rescaling, however, demonstrates that the histor- 
ical data themselves are largely discrepant with abso- 
lute values reaching up to an order of ma g nitude lower 
than the Allan data (see iHuege fc Falckel l|2003j) for a 
detailed discussion). It is therefore still imperative to 
gather independent data with good abs olute calibration 
with new experiments such as LOPES l|Horneffer et all 
2003). For air showers of energies around 10 17 eV this 
should be well feasible as the predicted absolute lev- 
els of emission \E(R,u>)\ around a few fiY m _1 MHz -1 
at 55 MHz are well above the galactic noise limit of 
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-0.4/0.15/0.05 fiV ra.- 1 MHz " 1 for a 3a detection w ith 
1/10/100 LOPES antcnna(s) l|Huege fc Falckei l2003|) . In 
areas with high radio-frequency interference levels such as 
the site of the KASCADE array, the noise levels are a 
factor of a few higher. 

Further major improvement of our model will be 
reached by the interfacing of our MC simulation to th e 
air shower simulation code CORSIKA l|Heck et al.Lll998j) . 
The solid particle physics foundation of CORSIKA will 
resolve the uncertainties related to the somewhat crude 
shower parametrisations us ed so far. (Detailed MC calcu - 
lations of the RICE group (|Razzaaue et all I2002L |2004|) . 
e.g., reveiled that the compactness of the particle distri- 
butions for electromagnetic showers in ice — and thus the 
coherence of the radio emission at high frequencies — was 
originally underestimated in analytical parametrisations.) 



Furthermore, our code provides a solid basis for the 
inclusion of additiona l effects such as Askaryan-type 
I As karva Cerenkov radiation and an inter- 

facing to the MC air shower simulation code CORSIKA. 
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Appendix: Trajectory 

Consider first a simple magnetic field geometry 



10. Conclusions 

We have successfully advanced our modelling of radio 
emission from cosmic ray air showers with elaborate 
Monte Carlo simulations in the time-domain. Our MC 
code takes into account the important air shower char- 
acteristics such as lateral and longitudinal particle distri- 
butions, particle energy and track length distributions, a 
realistic magnetic field geometry and the evolution of the 
air shower as a whole. The calculation retains the full po- 
larisation information and does not employ any far-field 
approximations . 

We predict emission patterns, radial and spectral de- 
pendences for an exemplary 10 17 eV vertical air shower 
and find good agreement with our earlier theoretical works 
and the historical data available. 

A major result that could not be obtained by analytic 
calculations alone is that asymmetries introduced into the 
total field strength emission pattern by the magnetic field 
direction are washed out completely in the radiation from 
an integrated air shower. Statistics of total field strengths 
alone can therefore not establish the geomagnetic emis- 
sion mechanism. The clear polarisation dependence on the 
magnetic field direction, on the other hand, allows a di- 
rect test of the geomagnetic emission mechanism through 
polarisation-sensitive experiments such as LOPES. 

After having documented the implementation details 
and having demonstrated the correctness and robustness 
of our MC simulations, our code is now in a stage where 
we can explore the dependence of the radio emission on 
a number of parameters such as shower axis direction, 
primary particle energy, depth of shower maximum and 
the like. Consequently, this will be our next step. 

Once these dependences are established, measure- 
ments of radio emission from cosmic ray air showers can 
be related directly to the underlying characteristics of the 
observed air showers. Due to the regularity and robustness 
of the modelled emission patterns, even a sparse sampling 
of the radiation pattern with a limited number of antennas 
would probably suffice for such an analysis. 



B = B 



(13) 



The (unperturbed) trajectory of a charged particle in a 
homogeneous magnetic field is a helix. Aligning the helix 
along the z-axis, it can be written as 



(-R B COS [u)B{t-t )) 
+R B sm[u B (t-t Q )} 
v cos a (t — to) 



where 



Rb 



vsma 



LOB 



(14) 



(15) 



denotes the radius of the projected circular motion and 
qeB 



UJ B 



7TO e C 



(16) 



is the gyration frequency of the particle with charge q 
times the elementary charge unit e and velocity 



c 



1 - 



1 



(17) 



The pitch-angle a is given by the (constant) angle be- 
tween the direction of the particle velocity vector and the 
magnetic field vector. 

To derive the general form of this trajectory, we first 
rotate the coordinate system such that the i?-field points 
in the desired direction. Afterwards, we adjust the phase 
to such that the particle's initial velocity has the desired 
direction as specified by the initial velocity vector V. In 
the last step we shift the trajectory so that at t = it 
coincides with the desired starting position R. 

We want to transform the simple geometry field B to 
the desired geometry 



cos 9 sin tp 
B = B | sin 9 sin ip 
sin 9 



(18) 
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where ip 6 [0, 2ir[ and 9 6 [0, n] are the azimuth and zenith 
angles known from spherical coordinates. This transforma- 
tion is achieved by applying a rotation matrix 
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B = T>B (20) 
and, inversely, 
B = D Y B. 

Applying the same rotation to the trajectory i|14|) yields 

r(t) = D f(t). (22) 

To infer the phase t corresponding to a given initial ve- 
locity V, we rotate back V to the simple geometry, 



V = D _1 V. 



(23) 



The x- and y-component of V then directly determine to 
through the relation 



oJBto = =F arctan 




(24) 



with the upper sign for q > and the lower for q < and 
where one has to take into account the correct quadrant 
for the arctan operation. In terms of the components of 
V this yields 



to = =p arctan 



± [(V x cos tp + V y sin p) cos 9 — V z sin 9] 



V y cos p — V x sin y 
The last operation that has to be employed is a translation 



r abs (i) = Ro + r(t) 
of the trajectory such that 
r abs (t = 0) =R + r(t = 0) 
which yields 
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The resulting trajectory r a b s (i) is thus fully defined for a 
given set of parameters R and V. The time-dependence of 
particle velocity and acceleration are then easily derived 
as the time-derivatives of r a b s (t). 



